Interactions

library(haven)
library(tidyverse)
library(marginaleffects) 
library(jtools)

states <- read_dta("states.dta")
nes <- read_dta("nes.dta")

# Supreme Court approval using feeling thermometer

# Don't forget to think about think about distribution of these DVs
ggplot(nes, aes(x=ft_SCOTUS)) +
  geom_histogram(color="white")

# We're interested in party and knowledge; see how they're coded, which is key for 
# interpretations of interactions
nes |> group_by(partyid7) |> summarize(n=n()) |> mutate(pct=100*n/sum(n))

nes |> group_by(polknow) |> summarize(n=n()) |> mutate(pct=100*n/sum(n))

# Model without interactions - unconditional effects of party and knowledge
summary(lm(ft_SCOTUS ~ partyid7 + polknow + Female + as.factor(Race3), data=nes))

# Model with interactions - conditional effects of party and knowledge
# First, use "full data" that includes missing data (lm takes it out automatically)
h <- lm(ft_SCOTUS ~ partyid7 + polknow + partyid7*polknow + as.factor(Female) + 
          as.factor(Race3), data=nes)
summary(h, digits=3)

# Alternative: If you include the multiplicative term, lm automatically
# includes the constituent terms. 
summary(lm(ft_SCOTUS ~ partyid7*polknow + as.factor(Female) + 
          as.factor(Race3), data=nes))

# Generate marginal effects of party id conditional on knowledge 
# We could do this manually -- write out equation

### Graph marginal effects (Brambor, Clark and Golder graph) 
# Use slopes command from marginaleffects package
mepid <- slopes(
  h, #model object
  variables = "partyid7", #variable for which we want marginal effect
  df = insight::get_df(h), #sets df for model to use t dist instead of z
  numderiv = "richardson", #SEs, delta method
  newdata = datagrid(polknow = 0:5)) #polknow the moderator, set range

# Use ggplot, dotted lines for 95% CIs; add horizontal zero line
ggplot(mepid, aes(x = polknow)) + 
  geom_line(aes(y = estimate)) +
  geom_line(aes(y = conf.high), linetype = 2) +
  geom_line(aes(y = conf.low), linetype = 2) +
  geom_hline(yintercept = 0) +
  labs(title="Marginal effect of Party ID", x="Political knowledge", 
       y="Marginal effect of Party ID")

# Use ribbon instead of dotted lines for 95% CI
ggplot(mepid, aes(x=polknow, y=estimate, ymin=conf.low, ymax=conf.high)) + 
  geom_line(color="black") +
  geom_ribbon(alpha=.3, fill="dodgerblue") +
  geom_hline(yintercept = 0, linetype="longdash", linewidth=.3) +
  labs(title="Marginal effect of Party ID", x="Political knowledge", 
       y="Marginal effect of Party ID")

# Use same concept, but use dots and rspikes
ggplot(mepid, aes(x=polknow, y=estimate, ymin=conf.low, ymax=conf.high)) + 
  geom_point(color="black", size=2) +
  geom_errorbar(width=0) +
  theme(plot.title = element_text(face="bold", size=12, hjust = 0.5)) +
  geom_hline(yintercept = 0, linetype="dashed") +
  theme_minimal() +
  labs(title="Marginal Effect of Party ID", x="Political knowledge", 
       y="Marginal effect of Party ID") +
  theme(plot.title = element_text(face="bold", size=12, hjust = 0.5))

# Flip to horizontal; typically used when you label the vertical axes with long labels
ggplot(mepid, aes(x=polknow, y=estimate, ymin=conf.low, ymax=conf.high)) + 
  geom_point(color="black", size=2) +
  geom_errorbar(width=0) +
  geom_hline(yintercept = 0, linetype="dashed") +
  theme_minimal() +
  coord_flip() +
  scale_x_continuous(breaks=c(0:5), 
                     labels=c('Very low knowledge','Low knowledge','Med-low knowledge',
                              'Med-high knowledge','High knowledge','Very high knowledge')) +
  labs(title="Marginal Effect of Party ID", x=NULL, 
       y="Marginal effect of Party ID") +
  theme(plot.title = element_text(face="bold", size=12, hjust = 0.5))

## Second way of graphing interactions: 
# Plot partial slopes for the effect of party on SC evals conditional on 
# low and high knowledge.

lknow <- predictions(
  h,
  by="partyid7",
  df = insight::get_df(h),
  numderiv = "richardson",
  newdata = datagrid(
    polknow = 0, partyid7=1:7, grid_type="counterfactual"))

hknow <- predictions(
  h,
  by="partyid7",
  numderiv = "richardson",
  df = insight::get_df(h),
  newdata = datagrid(
    polknow = 5, partyid7=1:7, grid_type="counterfactual"))

ggplot(data=lknow, aes(x=partyid7, y=estimate)) + 
  geom_line() +
  geom_line(data=hknow, aes(x=partyid7, y=estimate), color="dodgerblue") +
  geom_text(x=4, y=62, label="High pol. knowledge", color="dodgerblue") + 
  geom_text(x=5, y=53, label="Low pol. knowledge", color="black") +
  scale_x_continuous(breaks=c(1:7), labels=c('SD','WD','ID','I','IR','WR','SR')) +
  scale_y_continuous(limits = c(50, 65), breaks = seq(50, 65, by = 5)) + 
  theme_minimal() +
  labs(title="Conditional Effect of Party ID", x="Party ID",  
       y="Supreme Court Thermometer") +
  theme(plot.title = element_text(face="bold", size=12, hjust = 0.5))